source('00_util_scripts/mod_bplot.R')
source('00_util_scripts/mod_seurat.R')
library(tidyplots)
library(data.table)

# Gao 2021 JCI nasal wash, seq well ---------
jci.viro <-
  read_delim('trpm7_mucosa/data/GSE176269_CovidStudy_phenotype_061721.txt.gz')

jci.count <- Read10X_h5('trpm7_mucosa/data/GSE176269_CovidStudy_rawCounts.h5')

sobj <- jci.count |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj <- jci.viro |>
  mutate(.cell = ...1, .keep = 'unused') |>
  left_join(x = sobj, y = _)

sobj$orig.ident <- NULL

Idents(sobj) <- sobj$sampID

sobj$mito.ratio <- sobj |> PercentageFeatureSet('^MT-')

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj <- sobj |> NormalizeData()

sobj |> write_rds('trpm7_mucosa/gao21jci.rds')

sobj <- read_rds('trpm7_mucosa/gao21jci.rds')

sobj <- sobj |> filter(mito.ratio < 10)

sobj |> DotPlot('TRPM7', group.by = 'group', cols = 'RdYlBu')

sobj |>
  filter(group != 'unk_COVID') |> 
  DotPlot2d('TRPM7', status, group) |>
  pluck('data') |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM7 in human nasal wash',
       subtitle = 'Gao 2021 JCI')

sobj$status |> unique()

## M7 expr & IgG response --------
sobj |>
  filter(group == 'epithelial', cellType != 'unk_epi') |>
  DotPlot2d('TRPM7', status, cellType) |>
  pluck('data') |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', subtitle = 'Gao 2021 JCI',
       title = 'TRPM7 expression in epithelial cells from nasal wash')

sobj |>
  filter(group == 'epithelial') |>
  DotPlot('TRPM7', group.by = 'sampID', cols = 'RdYlBu')

epi.m7.exp <- last_plot() |>
  pluck('data')

igh <- sobj |> rownames() |> str_subset('IGH(A|D|E|M|G).{0,1}$')

sobj |>
  filter(group == 'Bcell') |>
  DotPlot(igh, group.by = 'sampID', cols = 'RdYlBu')

bc.igh.exp <- last_plot() |>
  pluck('data')

m7.order <- epi.m7.exp |>
  mutate(m7.expr = avg.exp.scaled, id, .keep = 'none')

epi.m7.exp <- epi.m7.exp |>
  bind_rows(bc.igh.exp) |>
  as_tibble() |>
  left_join(m7.order) |>
  mutate(id = fct_reorder(id, m7.expr, .na_rm = TRUE))

epi.m7.exp |>
  write_csv('epim7.b-igh.expr.covid.csv')

epi.m7.exp <- read_csv('epim7.b-igh.expr.covid.csv')

## IGHA1 & TRPM7 ------
epi.m7.exp |>
  filter(features.plot %in% c('IGHA2', 'TRPM7'), str_detect(id, 'Cov')) |>
  mutate(features.plot = ifelse(features.plot == 'IGHA2',
                                'IGHA2 (B cells)', 'TRPM7 (Epithelial cells)')) |>
  ggplot(aes(features.plot, avg.exp.scaled, group = id)) +
  geom_line() +
  geom_point(color = 'black') +
  labs(title = 'Average expression in COVID-19 nasal wash cells',
       y = 'Z-score scaled expression', x = 'Gene (Cell type)') +
  theme_pubr()

## IGHA1 & IGHA2 --------
epi.m7.exp |>
  filter(features.plot %in% c('IGHA1', 'IGHA2'), str_detect(id, 'Cov')) |>
  ggplot(aes(features.plot, avg.exp, group = id)) +
  geom_line() +
  geom_point(color = 'black') +
  stat_summary(aes(group = features.plot), geom = 'crossbar', fun = 'mean',
               width = .3) +
  labs(title = 'Average expression in COVID-19 nasal wash B cells',
       y = 'Normalized expression', x = 'Gene') +
  theme_pubr() +
  stat_compare_means(comparisons = list(c('IGHA1','IGHA2')), method = 't.test')

vlna12 <- sobj |>
  filter(group == 'Bcell', status == 'COVID-19') |>
  get_abundance_sc_long(features = c('IGHA1','IGHA2'))

vlna12 |>
  ggplot(aes(.feature, .abundance_RNA)) +
  geom_violin(aes(fill = .feature), scale = 'width') +
  stat_summary(fun = logtpm.mean, geom = 'crossbar',
               width = .3, color = 'black') +
  theme_pubr() +
  labs(x = 'Gene', fill = 'Gene', y = 'Normalized expression',
       title = 'Expression in COVID-19 nasal wash B cells') +
  stat_compare_means(comparisons = list(c('IGHA1','IGHA2')))

### covid ---------
epi.m7.exp |>
  filter(str_detect(id, 'Cov'), id != 'Cov10N') |>
  mutate(avg.exp.scaled = scale(avg.exp)[,1], .by = features.plot) |>
  BubblePlot() +
  labs(y = 'Individual', subtitle = 'Gao 2021 JCI',
       title = 'COVID-19 patients epithelial TRPM7 expression & B cell Ig expression')

m7hvl.cov.b.deg <- sobj |>
  filter(group == 'Bcell', status == 'COVID-19', sampID != 'Cov10N') |>
  FindMarkers(group.by = 'sampID', ident.1 = c('Cov14','Cov05', 'Cov15')) |>
  as_tibble(rownames = 'gene')

m7hvl.cov.b.deg |>
  plot_bill_volc(group1 = 'TRPM7-high', group2 = 'TRPM7-low')

igh <- igh |> str_subset('IGH[^E]')

sobj |>
  filter(group == 'Bcell', status == 'COVID-19', sampID != 'Cov10N') |>
  mutate(m7.type = sampID %in% c('Cov14','Cov05', 'Cov15')) |>
  bill.violin(igh, m7.type)

g1 <- last_plot()

g1 + theme(axis.text.x = element_blank()) +
  scale_fill_hue(name = 'Group',
                 labels = c('TRPM7-low (n=4)', 'TRPM7-high (n=3)')) +
  labs(x = 'Group', y = 'Normalized expression',
       title = 'B cell Ig expression in COVID-19 patients')

# flu -----------
epi.m7.exp |>
  filter(str_detect(id, 'f')) |>
  mutate(avg.exp.scaled = scale(avg.exp)[,1], .by = features.plot) |>
  BubblePlot() +
  labs(y = 'Individual', subtitle = 'Gao 2021 JCI',
       title = 'Flu patients epithelial TRPM7 expression & B cell Ig expression')

jci.viro |>
  filter(status == 'fluA') |>
  distinct(sampID, cDate, donorID)
